


sim_die_rolls <- function() {
  rolls <- sample(1:3,
                  size = 10,
                  replace = TRUE,
                  prob = c(1, 2, 3) / 6)
  
  c(sum(rolls == 1), sum(rolls == 2), sum(rolls == 3))
}



replicate_simulations <- function(n) {
  replicate(n, sim_die_rolls())
}


ncount <- 5000

results <- data.frame(t(replicate_simulations(ncount)))
colnames(results) <- c("one", "two", "three")

tm <- results %>%
  summarize(P=sum(one + two < 5) / ncount)


results %>%
  filter(two == 3) %>%
  summarize(one_m = mean(one))


dbinom(x=3, size=7, prob=1/4)



